Ground-state properties of the disordered Hubbard model in two dimensions 
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We study the interplay between electron correlation and disorder in the two-dimensional Hubbard 
model at half-filling by means of a variational wave function that can interpolate between Anderson 
and Mott insulators. We give a detailed description of our improved variational state and explain how 
the physics of the Anderson-Mott transition can be inferred from equal-time correlations functions, 
which can be easily computed within the variational Monte Carlo scheme. The ground-state phase 
diagram is worked out in both the paramagnetic and the magnetic sector. Whereas in the former a 
direct second-order Anderson-Mott transition is obtained, when magnetism is allowed variationally, 
we find evidence for the formation of local magnetic moments that order before the Mott transition. 
Although the localization length increases before the Mott transition, we have no evidence for the 
stabilization of a true metallic phase. The effect of a frustrating next-nearest-neighbor hopping t' 
is also studied in some detail. In particular, we show that t has two primary effects. The first 
one is the narrowing of the stability region of the magnetic Anderson insulator, also leading to a 
first-order magnetic transition. The second and most important effect of a frustrating hopping term 
is the development of a "glassy" phase at strong couplings, where many paramagnetic states, with 
disordered local moments, may be stabilized. 

PACS numbers: 71.30.+h, 71.27. +a, 71.55.Jv 



I. INTRODUCTION 

Within independent-electron approaches, all single- 
particle states are delocalized and the metallic or in- 
sulating behavior is determined by the existence of an 
energy gap between the highest occupied level and the 
lowest unoccupied one. However, whenever the Coulomb 
interaction becomes dominant over the kinetic energy, 
the independent-electron picture fails and electrons in 
the narrow bands close to the Fermi energy become lo- 
calized. Systems, whose insulating character is induced 
by electron correlations, are called Mott insulators^ The 
presence of disorder weakens the constructive interfer- 
ence that allows a wave-packet to propagate coherently 
in a periodic potential and may eventually lead to single- 
particle localization at the Fermi energy, hence to a fur- 
ther class of insulating materials, called Anderson insu- 
lators^ In this context, the conventional one-parameter 
scaling theory of conductance 3 - predicts that all solutions 
of the single-particle Schrcedinger equation in a disor- 
dered potential are localized in two dimensions (2D). 
Therefore, any amount of disorder in 2D drives a non- 
interacting electron system into an Anderson insulator. 
The inclusion of electron-electron interaction in weak 
coupling does not modify qualitatively the above conclu- 
sion^ although it may crucially affect physical properties 
like the tunneling density of states4£ Hence, it is widely 
accepted that 2D electron systems should always display 
an insulating behavior at sufficiently low temperature or 
large size no matter how weak the disorder is. 

Nevertheless, from time to time some indications have 
appeared that this conclusion might not be always cor- 
rect. Finkel'stein£ and Castellani et alZ considered the 
interplay between disorder and interaction by perturba- 



tive renormalization group methods and showed that, for 
weak disorder and sufficiently strong interactions, a 2D 
system might scale towards a phase with finite conduc- 
tivity. However, this result was not conclusive since the 
"metallic" region is found to occur outside the theory's 
limits of validity. Indeed, the common feature of these 
calculations is the crucial role played by spin fluctuations 
that grow large as the renormalization procedure is iter- 
ated. This tendency has been commonly interpreted as 
the emergence of local moments concomitantly with the 
progressive Anderson localization and, in continuum sys- 
tems, it might signal an incipient ferromagnetic instabil- 
ity^ In lattice models, this is likely to be substituted by 
a magnetic instability at some wave-vector determined 
by the topology of the Fermi surface. In this respect, 
an alternative approach could be to assume from the be- 
ginning a strong interaction that drives the system into 
a local moment regime, e.g., close to a magnetic Mott 
transition, and then turn on disorder. However, since 
it is already difficult to obtain an accurate description 
of a clean Mott transition, this approach is very hard 
to pursue. Moreover, assuming the scenario provided by 
dynamical mean-field theory (DMFT)^ one would im- 
mediately face with the so called Harris criterion^ In- 
deed, the correlation-length exponent predicted within 
DMFT for a clean Mott transition v ~ 1/2 is smaller 
than 2/D = 1 in D = 2 dimensions! 11 ' 12 This fact im- 
plies that the whole critical behavior has to be profoundly 
altered by disorder, which therefore can not be regarded 
at all as a weak perturbation. 

Only with the groundbreaking experimental work by 
Kravchenko and coworkers , 13 ' 14 the statement that lo- 
calization always occurs in 2D was really put into ques- 
tion. In fact, Kravchenko has been the first to observe 
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and claim that, above some critical carrier density, high- 
mobility silicon MOSFETs display a metallic behavior, 
i.e., a resistivity that decreases with decreasing tempera- 
ture down to the lowest accessible one. Below this critical 
density, the behavior of the resistance looks insulating, 
thus suggesting that a metal-insulator transition does oc- 
cur by varying the density. This experimental finding 
has renewed the interest in the interplay between disor- 
der and interaction and has generated further theoretical 
works along the same direction originally put forward 
by Finkel'stein, in the attempt to clarify some open is- 
sues and the applicability of the approachi 15 ' 16 However, 
the issue of having a genuine metallic phases in 2D is 
controversial and remains under debate, from both the 
experimental and the theoretical point of view.— ~— 

An approach where neither the interaction nor the dis- 
order are treated as a perturbation is therefore required 
to understand the complex physics of correlated disor- 
dered systems. The simplest model which contains both 
ingredients on the same level is the disordered Hubbard 
model, which has been intensively studied by several nu- 
merical methods, such as Hartree-Fock calculations in 
twc3 and three dimensions'^ extended DMFT'^ 
and quantum Monte Carlo simulations i^2r— However, 
all these methods have some drawback and only a com- 
bined analysis of complementary techniques can be able 
to clarify the nature of this challenging problem. Any 
approach based on single-particle descriptions, like un- 
restricted Hartree-Fock j 20 i 33 can uncover the emergence 
of an insulating gap only by forcing magnetic long-range 
order. More sophisticated approaches, like those based 
on DMFT ) 23 i 24 can in principle manage without mag- 
net ism ,25,26,34,35 but they usually miss important spatial 
correlations. The spatial distribution of local moments 
and its connection with Griffiths singularities, which 
may emerge close to the Mott insulator^ have been re- 
cently discussed within a Brinkmann-Rice approach of 
the Gutzwiller wave function.— 

In Ref. l38l an improved approach based on a varia- 
tional wave function has been proposed to deal simulta- 
neously with the physics of Anderson's and Mott's local- 
ization. In this work, we present an extensive and more 
complete study of the ground-state properties of the dis- 
ordered Hubbard model in two dimensions at half filling. 
Moreover, we also discuss in great details the variational 
method used in our calculations. Within our approach it 
is possible to describe both a direct transition between a 
compressible Anderson insulator and an incompressible 
paramagnetic Mott phase and a transition to a magnet- 
ically ordered insulator. In the first part, we discuss the 
variational wave function and we explain how it is possi- 
ble to distinguish between Anderson and Mott insulators, 
by means of static correlation functions. In the second 
part, we describe paramagnetic and magnetic transitions. 
By allowing for magnetic order, we show the evidence for 
the formation of local magnetic moments that order be- 
fore the Mott transition. Finally, we introduce a frustrat- 
ing (next-nearest-neighbor) hopping that favors a glassy 



behavior with many different local minima (and presum- 
ably long time scales). 

The paper is organized as follows: in Sec. [H] we present 
the method and the variational wave function that is used 
in our calculations; in Sec. IIIII we show the accuracy of 
our method; in Sec. lIVl we show our numerical results for 
the paramagnetic sector; in Sec.|Vl we present the results 
for the magnetic sector; finally, in Sec. I VII we draw our 
conclusions. 



II. MODEL AND METHODS 

A. Variational wave function 

We consider the two-dimensional Hubbard model with 
on-site disorder: 

^ = _ X! ^j C i,a C 3^ + h - C + '^2 e i n ^ + U ^,t n U W 
i,j,(T i i 

where c\ a (ci y(T ) creates (destroys) an electron at site i 
with spin a, and rij = J2 a n i," 1S ^ ne l° ca l density opera- 
tor, ej are random on-site energies chosen independently 
at each site i and uniformly distributed between [-D, D], 
U is the repulsive electron-electron interaction and t^ is 
the hopping amplitude. In the first part of this work, we 
consider only a nearest-neighbor hopping term t = 1 and, 
in the final part, we also add a frustrating next-nearest- 
neighbor term t' . We consider 45 degree rotated clusters 
with N = 2l 2 sites at half- filling, i.e., with N e = N elec- 
trons. 

For U = 0, the Hamiltonian (fl} reduces to the An- 
derson model for which the ground state is an Anderson 
insulator for any value of D, with gapless charge exci- 
tation but localized states. On the contrary, in the op- 
posite limit U/t — > oo, all charge fluctuations are sup- 
pressed, the system recovers translational invariance and 
the ground state becomes a Mott insulator. In this work, 
we study the zero-temperature properties for finite dis- 
order D and Coulomb interaction U , by using the vari- 
ational Monte Carlo algorithm. Our variational ansatz 
for the ground state is given by 

\*)=JG\SD), (2) 

where \SD) is an uncorrelated Slater determinant that is 
the ground state of a mean-field Hamiltonian: 

Hmf = - X] *y c !,<r c .7> + h - c - + X! £h<r n h°-> ( 3 ) 

i t j,a i.a 

where = t for neighboring sites and ti^ a are varia- 
tional parameters; in addition the next-nearest-neighbor 
hopping tij = t' is also used as a variational parameter 
in case of a finite value of t' in the Hamiltonian. We con- 
sider both paramagnetic and magnetic properties. In the 
former case, we impose the wave function | , F) to be para- 
magnetic, by fixing the variational parameters e,-^ = e^. 
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On the contrary, in order to study magnetic properties, 
we consider a variational wave function that may break 
the spin-rotational symmetry, namely we allow the vari- 
ational parameters to be e^j, ^ e^-j-. 

A more general choice of the trial wave function 
could also contain additional parameters for all (nearest- 
neighbor) hopping terms ty (with no translational invari- 
ance). Calculations performed on small lattices showed 
that, despite a much larger computational effort, this 
ansatz leads only to slightly lower energies, without mod- 
ifying the ground-state properties. Therefore, this possi- 
bility will be not considered in the following. 

The Gutzwiller factor Q is defined by 



g = exp[-^ffin?], 



(4) 



static structure factor N a 



N q = 



/doj , . 



(7) 



We emphasize that N q can be easily calculated within 
the Monte Carlo method through an equal-time density- 
density correlation over the ground state 



N q = (n- q n q ) - N, 



disc 



(8) 



where (n_ g n q ) = ($>\n- q n q \'ii) indicates the quantum 
average and the overbar denotes the disorder average. 
Notice that, in order to have a correct definition of the 
dynamical structure factor N q (uj), we have subtracted 
the disconnected term 



and J is a long-range Jastrow term, defined by 



J = cxp 



2 ^ 

1,3 



j(m - l)(rij - f) 



(5) 



While the Gutzwiller factors have been defined with a 
different parameter gi for each site in order to describe 
the non-homogeneous character of the system, we only 



consider translational invariant i>. 



1,3 



v(\r. 



|). In the 



following, the Fourier transform of the Jastrow param- 
eters will be denoted by v q . All these parameters can 
be optimized in order to minimize the variational energy. 
The choice of taking translational invariant Vij is done 
in order to reduce the total number of parameters and 
make the problem tractable from a numerical point of 
view. Moreover, since the Jastrow factor plays a primary 
role in the strong-coupling regime, where the disorder ef- 
fects are suppressed, this choice should not put serious 
limitations to our variational wave function. 

Summarizing, the variational parameters are i) the on- 
site energies ei j(T and, for t' ^ 0, the hopping t' of the 
mean- field Hamiltonian ([3]), ii) the Gutzwiller parame- 
ters gi, and iii) the translational invariant Jastrow pa- 
rameters Vij. The optimization of the variational state, 
without assuming any particular parametric form, is done 
by an energy minimization, which may deal with a large 
number (few hundreds) of parameters^ 



B. Correlation functions 

To work out the zero-temperature phase diagram, we 
make use of the f-sum rule that allows us to interpret the 
small-g behavior of density-density correlations. Starting 
from the definition of the average energy of the excita- 
tions 



J dujujN q (uj) 
J duN q (u>) ' 



(6) 



where N q (uj) is the dynamical structure factor. By in- 
tegrating over frequencies this quantity, we obtain the 



N^ c = (n_ q ){n q ), 



(9) 



which is related to the elastic scattering of electrons and, 
in a disordered system, is finite for generic momenta 
Notice that, in a disordered system q is no longer a good 
quantum number but, nevertheless, the average over dif- 
ferent disorder configurations restores the translational 
invariance. This fact suggests that the density-density 
structure factor N q can be a meaningful quantity to as- 
sess physical properties. 

After a straightforward calculation, we arrive to a sim- 
ple expression of A q for small momenta: 



lim A„ ~ 

9^0 9 N„ 



(10) 



From this equation, we have the important result that 
N q ~ \q\ implies the existence of gapless charge modes, 
since A q — » for \q\ — > 0, while if N q ~ q 2 charge ex- 
citations are presumably gapped, since A,j ~ const for 
| q | — > 0. These criteria have been already applied with 
success to the fermionic and bosonic Hubbard models 
without disorder; 41 ' 42 where it has been demonstrated 
that it is possible to describe a true Mott insulator by 
using a Jastrow-Slater wave function. In addition, it 
has been shown that there is a tight connection between 



N q and the Jastrow parameters 



For a clean sys- 



tem, a conducting state (metallic or superconducting in 
the fermionic case and supcrfluid in the bosonic one) is 
characterized by v q ~ I/|<?| in any spatial dimensions, 
whereas, in order to correctly reproduce the static struc- 
ture factor of a Mott insulator, v q must be more singular, 
i.e., v q ~ 1/q 2 in one and two dimensions and v q ~ 1/|<?| 3 
in three dimensions^ 

Let us now turn to the disordered model. Before con- 
sidering the interacting case, we would like to discuss 
the results for the non-interacting case and show that 
N q ~ \q\ is recovered, in agreement with the fact that 
the ground state is compressible. In Fig. [TJ we report 
the density-density correlations (n_ g n 9 ) and the discon- 
nected term Nq lsc calculated for a system with N = 800 
sites and disorder strength D/t = 5, averaged over 48 
disorder realizations. It is clear that both quantities are 
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FIG. 1: (Color online) Density-density correlation function 
{n- q nq} (upper panel), disconnected term N q ' BC (middle 
panel), and N q = {n- q n q ) — N q 1BC (lower panel) for the non- 
interacting Anderson insulator. The inset shows N q /\q\, from 
which it is clear that N q ~ |g|. Results are averaged over 48 
disorder realizations with D/t = 5 for N = 800 sites. 



finite for q — >• 0. However, once we consider the con- 
nected part of the density-density correlations N q , we 
have that N q ~ |g| for q — » 0, in agreement with the fact 
that the Anderson insulator is gapless, see Fig. Q] We 
would like to emphasize that, in contrast to the clean 
case where N q ~ \q\ implies a conducting behavior, here 
it just indicates a compressible state with gapless exci- 
tations, but not a metallic character, because the single- 
particle states are localized. 



III. ACCURACY OF THE WAVE FUNCTION 

The variational energy landscape of the disordered 
Hubbard model may be characterized by the presence 
of different local minima. In fact, if we start from differ- 
ent points in the parameter space, namely from different 
values of g^, Vij and 6j j0 -, we may converge to different so- 
lutions. In the following, whenever we consider the para- 
magnetic sector, we impose — along the whole 
optimization procedure. On the other hand, when mag- 
netic solutions are allowed, these conditions are relaxed 
along the Monte Carlo simulation and the two on-site 
energies (for up and down spins) are optimized indepen- 
dently. In the latter case, the starting configuration can 
be taken to be either paramagnetic (i.e., with Sj iT = e^j.) 
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FIG. 2: On-site magnetization (m;) for a typical disorder con- 
figuration with D/t = 5, U/t = 10 (left panels) and U/t = 16 
(right panels). The upper panels correspond to the wave func- 
tion obtained starting from a paramagnetic point, whereas the 
lower panels correspond to the solution obtained from a stag- 
gered point. The black contour shows the elementary cell of 
the lattice which is repeated to mimic the infinite lattice with 
periodic boundary conditions. 
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FIG. 3: Accuracy of the variational energies AE = (Eo — E v ) 
(where Eo and E v are the exact and the variational energies, 
respectively) for different wave functions on a 4 x 4 lattice 
with D/t = 5: the correlated magnetic state (full squares), 
the paramagnetic state (empty circles), the magnetic state 
without Jastrow factors (empty squares), the Hartree-Fock 
state (full circles). The exact ground-state energy is computed 
by the Lanczos algorithm. 



or with a small staggering of the magnetization (i.e., with 
(i.a = Ci+a(— l)\ Xi+Vi \6, where S is a small quantity). We 
would like to stress that, even by considering a paramag- 
netic starting point, the converged solution will generally 
have e i)T ^ e^. 

In contrast to the paramagnetic case, in which the en- 
ergy landscape usually has one minimum (i.e., the same 
parameters are obtained when starting from different ini- 
tializations), when allowing a magnetic wave function dif- 
ferent local minima may appear. This feature is par- 
ticularly evident for large enough Coulomb repulsion, 
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whereas in the weak-coupling regime we recover a sim- 
ple picture with only one minimum. Remarkably, the 
appearance of different local minima is related to the 
presence of short-range magnetic correlations. In fact, 
by increasing the on-site Coulomb repulsion, some sites 
acquire a finite magnetization and eventually order, giv- 
ing rise to the typical staggered pattern. Whenever lo- 
cal moments are present or the magnetization is very 
small, there are different electronic arrangements that 
give similar energies but may be hardly connected by 
simple single-particle moves, so to define metastable lo- 
cal minima. However, especially for the unfrustrated or 
the weakly-frustrated case, the presence of these local 
minima is not a dramatic problem. In fact, generally, 
all physical quantities, as for instance the density-density 
structure factor N q , give similar results in all these cases. 
Therefore, we can safely conclude that all the optimized 
states share the same physical properties. As an example, 
Fig. [2]shows the on-site magnetization (mj) = {rii f— ^4,) 
pattern for the variational wave function optimized both 
starting from the paramagnetic point and from the stag- 
gered point. It is evident that there is no considerable 
difference between the two states, although some of the 
magnetization values are slightly different. 

We remark that, in most cases, we obtain a lower vari- 
ational energy by starting from the paramagnetic point, 
even if the final converged state is magnetically ordered. 
Therefore, although in some particularly delicate cases 
we considered different starting points, we usually ini- 
tialize the simulation with a paramagnetic configuration. 

Let us now discuss the accuracy in energy for a 4 x 4 
lattice, where the exact ground state can be calculated 
exactly by the Lanczos algorithm. In particular, we con- 
sider four different wave functions: i) the magnetic state 
with on-site Gutzwiller and Jastrow factors (that corre- 
sponds to our best ansatz), ii) the paramagnetic state 
with Gutzwiller and Jastrow terms, iii) the magnetic 
state with only Gutzwiller projectors, and iv) the mag- 
netic mean- field state \SD) (i.e., without any correla- 
tion term). For U = 0, the exact ground state wave 
function can be obtained in all these cases, implying a 
very good accuracy also for small but finite values of 
U/t. For small interactions there is no appreciable dif- 
ferences between paramagnetic and magnetic wave func- 
tions and for all correlated states the accuracy in the 
energy, i.e., (E — E v )/E (where E and E v are the ex- 
act and the variational energies, respectively), is lower 
than 1%. However, even for U/t = 4, the Hartree- 
Fock state, with no Gutzwiller and Jastrow factors, give 
a much worse accuracy than the other three correlated 
wave functions, see Fig. [3] For larger values of the in- 
teraction U, the situation is different, since the param- 
agnetic state is generally worse than the magnetic wave 
functions, including the Hartree-Fock state. The accu- 
racy of our best ansatz is about 10% (or less) up to very 
large Coulomb repulsions, which is acceptable in a dis- 
ordered model. However, we notice that the long-range 
Jastrow factor is not crucial, and wave functions i) and 



iii) give comparable energies both for small and large in- 
teraction values, e.g., U/t — 4 and U/t = 16, see Fig. [3] 
In fact, on the one hand, the on-site Gutzwiller factor can 
easily account for the small charge correlations induced 
by the Coulomb repulsion in the weak-coupling regime. 
On the other hand, for large U/t, the ground state has 
strong magnetic correlations and is well described with 
a (gapped) mean-field state. In the intermediate regime, 
our magnetic state with both on-site Gutzwiller projec- 
tors and the long-range Jastrow term may give a con- 
siderable improvement over the other states considered 
here. It should be stressed, however, that the Jastrow 
factor is essential to have a paramagnetic Mott insula- 
tor, since in this case the charge gap cannot be opened 
by an uncorrelated Hartree-Fock state. On the contrary, 
in the magnetic case, long-range Jastrow correlations are 
not strictly necessary to capture the correct nature of the 
ground state, since the charge gap can be naturally cre- 
ated by a finite antiferromagnetic mean-field parameter. 



IV. RESULTS: THE PARAMAGNETIC CASE 

Let us start our analysis of the disordered Hubbard 
Hamiltonian (jTJ) by enforcing a paramagnetic variational 
wave function, e.g., by fixing the constraint e^-j- = ej j. 
for the parameters of the auxiliary mean-field Hamilto- 
nian ([3]). Although this choice is biased, since antiferro- 
magnetic order may be present for finite electron-electron 
interaction, it gives significant insights into the interplay 
between disorder and interaction, without any "spuri- 
ous" effect due to magnetism. We apply the f-sum rule 
to distinguish the compressible Anderson insulator from 
the incompressible Mott insulator, i.e., we look at the 
different behavior of N q and v q for different values of 
interaction U and disorder D. Here, we take a rather 
strong disorder (D/t = 4, 5, and 6) in order to have a 
localization length that is smaller than the numerically 
accessible system sizes. 

A detailed analysis of these quantities allows us to 
identify the Mott transition at C/ C MI = (11.5 ± 0.5) t 
for D/t = 5M- Indeed, for small values of the interac- 
tion strength, i.e., for U < U^ 11 , we have that N q ~ 
whereas N q ~ q 2 in the strong-coupling regime U > U^ 11 . 
The latter behavior is symptomatic of the presence of 
a charge gap hence of a Mott insulating behavior4i We 
would like to emphasize that the present results are qual- 
itatively similar to those obtained within the clean Hub- 
bard model, with linear coefficient of N q going smoothly 
to zero at the phase transition, indicating that the tran- 
sition is likely to be continuous. We mention that the 
Fourier transform v q of the optimized Jastrow parame- 
ters Vij is also compatible with a Mott transition at the 
same value of the interaction strength, with its small-q 
behavior changing from v q ~ l/\q\ to v q ~ I/9 2 across 
the Mott transition^ 

Most interestingly we find that iV^ lsc goes continuously 
to zero at the phase transition, see Fig. 2] This fact 
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FIG. 4: (Color online) Disconnected part of the density- 
density correlation function Nq' BC for the smallest q as a func- 
tion of the interaction U for D/t = 5 and N = 98 sites. 
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FIG. 5: (Color online) The same as in Fig. [4] for different 
momenta. 



allows us to identify a simple and variationally accessible 
quantity to distinguish between an Anderson insulator 
and a Mott insulator. We would like to stress that this 
behavior is not restricted to q — > 0, but we recover a 
similar trend also for finite momenta, although disorder 
fluctuations are larger for larger q vectors, see Fig. [5] 
These results demonstrate that charge fluctuations are 
strongly suppressed and eventually become local as U/t 
increases. 

All these results suggest that disorder is strongly sup- 
pressed in the regime of strong correlations. This fact is 
also corroborated by the calculation of the variance of the 
distribution of the on-site energies ej i<T , which was shown 
in Ref.fH. Previous Monte Carlo calculations have shown 
that a repulsive interaction may screen the local energies, 
thereby generating an effectively weaker random poten- 
tial^ At the mean-field level^ and beyond^ the screen- 
ing effects have been widely discussed. The Hartree-Fock 
state leads to a disorder screening only for moderate in- 
teractions, while it gives almost unscreened on-site ener- 
gies in the strongly-correlated regime. On the contrary, 
our correlated variational approach is able to capture the 
correct physics also for large interaction U, where the 
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FIG. 6: (Color online) Static structure factor N q divided by 
\q\ as a function of |q|, for different values of the interaction 
U. Upper panel: U/t — 2 (circles), 3 (upward triangles), 4 
(squares), and 5 (downward triangles). Lower panel: U/t = 
6 (circles), 7 (upward triangles), 8 (squares), 9 (downward 
triangles). Calculations have been done for D/t — 5 and 
N = 98 sites. 
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FIG. 7: Phase diagram for the disordered Hubbard model in 
the paramagnetic sector. 



disorder, though finite, is highly suppressed. The redis- 
tribution of on-site energies leads to a decreased localiza- 
tion of the electronic state at the Fermi level. Neverthe- 
less, the single-particle eigenstates are always localized, 
even though a very large localization length may develop. 
Within our variational approach, a important question is 
whether the action of the Gutzwiller correlator and the 
Jastrow factor could turn a localized \SD) into a delocal- 
ized many-body state |\&). Unfortunately, the variational 
method does not give access to dynamical quantities and, 
therefore, we cannot make a definite statement. Never- 
theless, we tend to believe that such a transmutation of 
a localized \SD) into a delocalized is unlikely. In any 
case, the previous results show an increase in the "metal- 
licity" of the ground state, with a partial screening of 
disorder. This result is in agreement with the fact that 
the linear slope of N„ has a non-monotonic behavior as a 



7 



12 
I 8 
a. 4 



U/t=3 









40 










U/t=10 — 


32 




U/t=11 — - 


24 






16 




1:::: 


8 



Li 





32 
^ 24 

w; 16 

Q_ 



U/t=6 



0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 
150 

100 
50 





0.0 0.1 0.2 0.3 0.4 0.00 0.05 0.10 0.15 

£|(q) 



£j(q) 



FIG. 8: (Color online) Distribution of £j(g) evaluated at q 
(n/7,n/7). Calculations are done for D/t = 5 and N — i 
sites. 
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FIG. 9: (Color online) The same as in Fig.[8]but for the local 
gap Aj(q). 



function of U, showing a peak for U/t ~ 7 that indicates 
an accumulation of low-energy states around the Fermi 
energy, see Fig. HI In fact, the linear slope of N q is re- 
lated to the compressibility of the system. Nevertheless, 
it has to be noticed that the Slater determinant \SD) is 
the ground state of a mean- field Hamiltonian that always 
describes non-interacting electrons with on-site disorder, 
no matter how large the Coulomb repulsion is. Therefore, 
even though the single-particle eigenstates may have a 
very long localization length because of the suppression 
of the effective on site disorder, yet this length remains 
finite in two dimensions. 

The analysis of the density-density correlation func- 
tion, the Jastrow factor, and N^ lsc for different values of 
D/t allows us to draw the paramagnetic phase diagram 
in the (U,D) plane, see Fig. [?l 

Finally, in order to gain a deeper understanding on 
the local behavior, namely how each single site behaves 



across the Anderson-Mott transition, we introduce a local 
f-sum rule: 



(11) 



Here Nj(q) defines the local static structure factor: 

N M = I? £<"* n J w^*"^' (12) 



where {ninj) conn represents the quantum average of riiTij 
after subtracting the disconnected term (m){nj}. In ad- 
dition, Sj (q) is related to the local kinetic energy 



-It 



E 



(c] Cji(T + h.c.)(e^- R ^-l) (13) 



where indicates the sum over the nearest neighbors of 
the site j. Both T,j(q) and Nj(q) can be easily evaluated 
in the variational Monte Carlo scheme, since they require 
the computation of equal-time correlations. 

The limit for small momenta of Eq. (jlip gives impor- 
tant insights into the local gap, making it possible to 
understand if at the Mott transition all sites become lo- 
calized simultaneously, or non-homogeneous fluctuations 
are still present. In Fig. [8j we report the distribution of 
Tij(q), evaluated at the smallest value of the q available 
within a 98-site lattice, i.e., q = (iy/7,tt/7). The average 
value of slightly increases by increasing U and it 

has a maximum in the regime U ~ D (similarly to what 
happens to N q /\q\, see Fig. [5]). The distribution P[Ej(q)] 
is rather large for small U/t and shrinks when the in- 
teraction strength is increased, in agreement with the 
fact that disorder is suppressed by interaction. However, 
even very close to the Mott transition, the variance of 
P[T,j(q)] stays quite large, indicating that a considerable 
number of sites still has large fluctuations. This fact can 
be interpreted as a two-fluid behavior, where a fraction 
of sites can be regarded as localized particles, whereas 
the remaining ones behave like in the Anderson insula- 
tor. Therefore, the Mott transition is driven by only a 
fraction of the total number of sites. For U > U^r , we 
recover a situation where all sites can be ascribed to the 
Mott phase and P[Ej(q)} has a very sharp peak, with a 
small variance. The distribution of Nj (q) is very similar 
to the one of the local kinetic term. By contrast, the 
distribution of the local gap Aj(g) is rather narrow for 
U < U^, where for small momenta Aj(q) ~ 0, due to 
a vanishing gap in the Anderson phase, see Fig. [9j Nev- 
ertheless, the distribution has very long tails (with very 
small weight), which are related to disorder fluctuations; 
these tails tend to be suppressed by increasing the inter- 
action U. For U > Uc , a charge gap opens up in the 
average density of states but the size of the gap turns out 
to be different from site to site, which implies a rather 
broad distribution, see Fig. [SJ 
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V. RESULTS: THE MAGNETIC CASE 

On the square lattice at half filling, in the absence of 
frustration and disorder, an arbitrarily weak repulsive 
Hubbard interaction U is able to induce long-range anti- 
ferromagnetic order. In fact, in this case, the presence of 
a perfect nesting in the Fermi surface implies a diverg- 
ing susceptibility at Q — (tt, it) that, in turn, opens a 
finite gap at the Fermi level. Therefore, the ground state 
is a band insulator for any finite value of the interaction 
U > 0. By contrast, in the presence of a local random po- 
tential, the charge gap may be filled by (localized) energy 
levels, possibly destroying the long-range magnetic order. 
Here, we address the important problem of the competi- 
tion between Anderson localization and magnetic order 
by using our improved variational approach and allowing 
for a magnetic Slater determinant \SD) with ^ i^. 
First, we discuss the phase diagram of the disordered 
Hubbard model ([1} with only a nearest-neighbor hop- 
ping t. In this case, a finite value of the interaction U^ F 
is needed to have the onset of long-range magnetic order: 
below U^ F the system is described by a standard param- 
agnetic (compressible) Anderson insulator, above U^ F a 
finite antiferromagnetic order parameter develops; how- 
ever, the excitation spectrum remains gapless and the 
system is compressible. By further increasing the inter- 
action U, i.e., for U > U^ 1 the ground state undergoes 
a second phase transition to an incompressible antiferro- 
magnetic insulator with a finite charge gap. Interestingly, 
in the paramagnetic Anderson insulator, local moments 
with a finite value of (m,) = {rii^ — n^) develop, sug- 
gesting that itinerant electrons may not be able to fully 
screen magnetic impurities created by disorder. 

In the last part, we add a next-nearest-neighbor (frus- 
trating) hopping t' . Also in this case, we show that the 
Mott insulating phase is always accompanied by mag- 
netic order, although with a sufficiently large ratio t'/t 
many local minima appear in the energy landscape, with 
competing magnetic properties. In particular, we find 
that the lowest-energy solution displays magnetic long- 
range order but many other disordered states with local- 
ized moments may be stabilized. 



A. Magnetic phase diagram 

Let us now consider large systems. In order to assess 
the magnetic properties, we define the total magnetiza- 
tion: 



(14) 



where m, = n^-j- — n^. In analogy with the clean model, 
also in presence of disorder, by increasing the electron- 
electron repulsion, there is a tendency toward magnetic 
order at Q = (tTjTt), and, therefore, we concentrate on 
this value of the momentum. 
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FIG. 10: (Color online) Staggered magnetization M and the 
disconnected term of the density-density correlations Nq lac as 
a function of U for disorder D/t — 5 (upper panel) and as a 
function of D for U/t = 4 (bottom panel). All calculations 
have been done for N = 98 sites. 
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FIG. 11: Phase diagram for the disordered Hubbard model in 
the magnetic sector. 



In Fig. [TOJ we report our results for D/t = 5 and dif- 
ferent values of the Coulomb repulsion and for U/t = 4 
and various disorder strength. Fixing D/t = 5, we find 
that U^ F = (6.5 ± 0.5) t, since the magnetization is fi- 
nite for U > C/ AF whereas it vanishes for U < U^ F . 
Moreover, as we discussed in the previous section, the 
information about charge fluctuations can be worked out 
from either N q or the disconnected term iV^ lsc . We ob- 
tain that iV 9 disc 



vanishes at £/ c MI = (10.5 ± 0.5) t, sig- 
naling the opening of a charge gap, see Fig. [TDJ Remark- 
ably, there is a finite region in which both the magnetiza- 
tion and the compressibility fluctuations are finite. This 
fact implies a stable regime that shows antiferromagnetic 
long-range order without a charge gap. We notice that 
this intermediate phase is much reduced when consider- 
ing U/t = 4 and vary the disorder strength, see Fig. [TO] 
In this case, we can estimate that D^ F = (2.5 ± 0.5) t 
and D™ 1 = (1.5 ± 0.5) t. These results lead to the phase 
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tained with the full variational wave function (full circles), 
containing both Gutzwiller and Jastrow terms, and the one 
obtained with no Jastrow factor but only the Gutzwiller pro- 
jector (full triangles). N q is shown in the upper panel and 



in the lower panel. 



diagram sketched in Fig. [TT] For U = the system is a 
(paramagnetic) Anderson insulator for every finite disor- 
der D > 0. Instead, for D = the ground state is a Mott 
insulator with antiferromagnetic order for every U > 0. 
When both disorder and interaction are finite, there is an 
intermediate phase between the paramagnetic Anderson 
insulator and the antiferromagnetic Mott insulator. This 
phase is characterized by long-range magnetic order, but 
also by a finite compressibility. Although some authors 
identified this phase with a metal^ we do not find any 
evidence in favor of it (see below). 

Let us now consider in more detail the nature of the 
Anderson-Mott transition that emerges from our varia- 
tional approach, once we allow for spin-rotational sym- 
metry breaking. We would like to remind the reader 
that, in the paramagnetic case, the Mott insulator can 
be obtained only thanks to a singular Jastrow factor, 
i.e., v q ~ 1/q 2 - In this case, the (paramagnetic) mean- 
field Hamiltonian ^ is always gapless and the charge gap 
opens because of the strong correlations induced by the 
Jastrow term. In the magnetic case instead, two different 
mechanisms can open a gap: i) the long-range charge cor- 
relations induced by the Jastrow factor and ii) the onset 
of long-range antiferromagnetic order due to a staggering 
of the e^'s. For U < £/ c MI = (10.5 ± 0.5) t, the static 
structure factor behaves like N q ~ \q\ and the Fourier 
transform of the Jastrow parameters like v q ~ l/|g|; on 
the other hand, for U > U^, we have that N q ~ q 2 and 
v q r~j 1/q 2 . Thus in the intermediate phase with long- 



range magnetic order and finite compressibility, N q ~ \q\ 
and v q ~ l/\q\- In order to understand which is the 
most relevant ingredient that opens the charge gap, we 
calculate the structure factor N q and the disconnected 
term N q lsc close to Mott transition for the full variational 
wave function |\I>) = JQ\SD) and for another state that 
only contains Gutzwiller terms, i.e., = Q\SD). The 
Slater determinant \SD) is independently optimized in 
the two cases. The results for N q and N q lsc are reported 
in Fig. [T2] and indicate that the behavior for the two 
states is very similar; even the critical value U^ 11 for the 
Mott transition is the same in the two cases. Therefore, 
we can conclude that, in the magnetic case, the charge 
gap opens mainly because of the presence of the mean- 
field order parameter. However, within the correlated 
wave function the Jastrow parameters still behave 
like v q ~ 1/q 2 in the Mott phase. 

In summary, the following scenario emerges: in the in- 
termediate phase with antiferromagnetic order but finite 
compressibility, the mean-field density of state is large 
at the Fermi level, however all states are localized; by 
increasing the interaction U, the Jastrow factor becomes 
stronger and, at the same time, there is a suppression 
of the mean-field density of states at the Fermi level; by 
further increasing U, a single-particle gap opens and the 
system becomes incompressible, see Fig. [T3J We remark 
that the tendency towards metallicity for intermediate 
values of U/t is suppressed by the presence of magnetic 
order: this can be seen by noticing a reduced density of 
states at the Fermi level for U/t ~ 8. In fact, for small 
values of the interaction, the localization length is short 
because of the strong disorder; then it becomes larger 
for higher values of the interaction U due to the disorder 
screening and then, at the antiferromagnetic transition, 
it decreases again. Moreover, we notice that the single- 
particle wave functions are more localized in the magnetic 
case than in the paramagnetic one, even in the regime of 
maximum "delocalization" , i.e., U ~ D. Unfortunately, 
a full analysis of size scaling of the localization length 
(by considering the inverse participation ratio, as used in 
Ref. [20h is very hard, because the available sizes do not 
allow us to reach definitive conclusions. 

Let us now focus on the development of local mag- 
netic moments in presence of the Coulomb interaction, 
in order to demonstrate that they may appear before the 
metal-insulator transition. For this issue, we consider the 
following quantity, which is related to the fluctuations of 
the local magnetization 



M L = 



(15) 



In a paramagnetic state with local moments, namely a 
state in which some sites have an on-site magnetization 
(m/) =/= 0, the total staggered magnetization is vanishing, 
i.e., M = 0, while Ml is finite. On the contrary, in the 
antiferromagnetic phase Ml ~ M. Therefore, by com- 
paring Ml and M, it is possible to have a good feeling 



10 




go 

Q 0.10 



FIG. 13: (Color online) Evolution of the density of states 
(DOS) of the auxiliary mean-field Hamiltonian as a function 
of the interaction U for N = 98 sites. The case with spin 
dependent on-site energies ?i )CT is considered. 
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FIG. 15: (Color online) Probability distribution of the of the 
absolute value of the on-site magnetization |(mi)| for N = 98, 
D/t — 5 and different values of the interaction U/t. 
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FIG. 14: (Color online) Staggered magnetization M (blue tri- 
angles) and fluctuations of the local magnetization Ml (red 
circles) as a function of U for disorder D/t — 5 (upper panel) 
and as a function of D for U/t = 4 (bottom panel). Calcula- 
tions have been done for N = 98 sites. 



on-site magnetization (mj) (as shown in Ref. 1 3 81) or from 
the probability distribution of |(mi)|, see Fig. [TS] Here, 
we observe that for small interaction values, U/t ~ 3, 
the distribution has a narrow peak in correspondence of 
|(mj)| = 0, however, at the same time, it has long tails 
indicating the presence of local moments. In this regime, 
the ground state is an Anderson insulator with a large 
number of paramagnetic sites and (n,-) = 0,1,2. For 
U ~ U^ F the distribution is spread between and 0.8, 
with a peak in correspondence of ~ 0. This fact 

highlights the coexistence of paramagnetic sites with lo- 
cal magnetic moments; these sites are not spatially cor- 
related hence long-range magnetism is absent. By in- 
creasing the interaction strength, the peak at ~ 
disappears and the one at |(77ij)| ~ 1 becomes more pro- 
nounced; in this case, local moments eventually display 
the typical staggered pattern of Neel order. Nevertheless, 
charge excitations are still gapless and N q ~ Finally, 
in the Mott insulating phase the distribution has a nar- 
row peak at |(/nj)| = 1. 



B. Frustrating case 



on the presence of local moments in the ground state. 
In Fig. [HI we show the staggered magnetization M and 
Ml both for D/t — 5 and different values of the interac- 
tion U and for U/t — 4 and different disorder strengths 
D. We find that for U > C/ C AF = (6.5 ± 0.5) t [and for 
D < D^ F = (2.5 ± 0.5) t] the two magnetization val- 
ues are very close, while in the paramagnetic phase we 
observe that Ml > M ~ 0. This fact suggests a magneti- 
cally disordered phase in which the on-site magnetization 
(mi) is finite for some sites. We identify those sites with 
(m;) ^ as local magnetic moments. The existence of 
such moments can be extracted from the pattern of the 



In the previous section, we showed that a disordered 
system of electrons on a square lattice undergoes a mag- 
netic transition before becoming a Mott insulator. There- 
fore, the Mott insulator is generically accompanied by 
magnetic order. However, in disordered materials, one 
could expect that long-range order may be strongly sup- 
pressed, leading to a bona-fi.de Mott transition, where 
the incompressible phase has no magnetic order. In this 
sense, the presence of a next-nearest-neighbor hopping 
t! may help to approach the Mott phase without any 
"spurious" magnetic effects. Indeed, this kind of frus- 
trated hopping generates, in the strong-coupling regime, 
a super-exchange term J' that competes with the nearest- 
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FIG. 16: (Color online) Collinear magnetization M (upper 
panel) and the disconnected term of the density-density cor- 
relations Nq lsc (lower panel) of the lowest energy solution as a 
function of U for D/t = 5 and t' jt = 1. The shaded area indi- 
cates the region where paramagnetic and magnetic solutions 
have similar energies. 



neighbor one J. 

First of all, we notice that a finite frustrating ratio 
t'/t generates very complicated energy landscapes, with 
many local minima. Furthermore, in contrast to the un- 
frustrated model, where different local minima share very 
similar physical properties, here different starting points 
in the parameter space may give rise to rather differ- 
ent wave functions, especially in the intermediate and 
strong-coupling regimes. Concerning the starting point 
for the energy optimization, we will consider i) a para- 
magnetic point, with e^-j- = e^j., ii) a staggered point, 
with €i t<J = £i + a(— ly Xi+ y*\d, and hi) a collinear point, 
with 6j jCT = rr(— ly x *\ . Notice that, in the latter case, 
the starting choice breaks both translational and rota- 
tional invariance, so to favor collinear magnetic order 
with Q = (n, 0), which is suitable for large t'/t. Simi- 
larly, we can also consider e ij(T = cr(— ly Vi ', which gives 
rise to a collinear order with Q — (0,tt). In all these 
cases, each local energy ej l<T (for up and down spins) is 
independently optimized, in order to achieve a full energy 
minimization. 

For weak and intermediate frustration, the outcome is 
rather similar to the case with t' = 0, except for U ~ 
U^ F . For instance, for t'/t = 0.6 and D/t = 5, the para- 
magnetic Anderson insulator is stable for U < U^ F ~ 6t. 
Up to this value of the Coulomb interaction, the energy 
of the converged state does not depend upon the choice 
of the starting point (as for the unfrustrated case) and 
there is no evidence for any non-trivial magnetic pattern. 
On the contrary, for 6 < U/t < 12, different results for 
the staggered magnetization are found when we initial- 
ize according to i) or ii). In particular, the magnetiza- 
tion is strongly suppressed when considering a param- 
agnetic initial condition. These different solutions are 



very close in energy up to U ~ 8t, even though, usually, 
the lowest minimum shows long-range order and finite 
compressibility; for larger values of the interactions, i.e., 
8 < U/t < 12, the magnetic state gives definitely the best 
energy. By further increasing the Coulomb repulsion, 
namely for U > U^ 11 ~ 12t, the system becomes a Mott 
insulator with Neel order at Q = (tt,tt). These results 
suggest that, for a finite frustrating ratio, the magnetic 
transition may become weakly first order (with a regime 
in which there is a coexistence of different phases) . 

The existence of a finite region with nearly degenerate 
states with different magnetic properties is even more 
pronounced in the strong frustration regime. For t' = t, 
the solution obtained starting from ii) gives always an 
higher energy with respect to i) and iii), demonstrating 
that the Neel state is clearly disadvantaged. Again, in the 
weakly-correlated regime, i.e., U < U^ F ~ 8 t, paramag- 
netic and magnetic starting points give similar energies 
and magnetization patterns. For U > U^ 11 ~ 13 t, the 
ground state is a gapped insulator with A^ lsc = 0, see 
Fig. [THl Instead, in the intermediate region, for U^ F < 
U < U^ 11 , the energy landscape shows different minima, 
which are nearly degenerate up to U ~ lit, whereas, for 
larger interactions, the wave function with collinear order 
gives the best approximation for the ground state, with 
a clear long-range magnetic order, see Fig. [T6l 

The remarkable feature is that, in a wide range of 
Coulomb repulsions, it is possible to find low-energy 
states just by starting from a paramagnetic wave func- 
tion. This choice gives rise to patterns in which most of 
the sites have a net magnetization but an overall vanish- 
ing magnetic order. For U/t ~ 16, these solutions are 
incompressible, i.e., N^ c ~ and, therefore, may be 
viewed as disordered Mott insulators. By decreasing the 
interaction strength, these states turn compressible, still 
having a large number of local moments. The presence of 
these metastable solutions, which are almost degenerate 
with the magnetically ordered wave function, suggests a 
sort of "spin-glass" behavior. The existence of a large 
number of such disordered states prevents one to have 
a smooth convergence to the lowest-energy solution, by 
starting from a generic configuration. 

In summary, the frustrating hopping t' has two pri- 
mary effects. The first one is the narrowing of the stabil- 
ity region of the magnetic Anderson insulator. In addi- 
tion, we have evidence that the magnetic transition turns 
to be first order, in contrast to the unfrustrated case. 
The second and most important effect of a frustrating 
hopping term is the development of a "glassy" phase at 
strong couplings, where many paramagnetic states, with 
disordered local moments, may be stabilized. Although 
we do not find any evidence in favor of the stabilization of 
a true (non-magnetic) Mott phase, this possibility cannot 
be ruled out in the whole phase diagram. 
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VI. CONCLUSIONS 

In this paper, we have studied, by means of a varia- 
tional Monte Carlo technique, zero-temperature proper- 
ties of the disordered two-dimensional Hubbard model 
at half filling. First of all, we showed that a varia- 
tional wave function is able to describe the Anderson- 
Mott transition without any symmetry breaking, i.e., a 
transition from a paramagnetic Anderson insulator to a 
paramagnetic Mott insulator. This is achieved thanks to 
long-range charge correlations induced in the wave func- 
tion by a Jastrow factor. We showed that the transition 
can be easily detected within variational Monte Carlo 
by looking at the behavior of the static structure factor 
and of the Fourier transform of the Jastrow parameters, 
namely following the same criteria for the Mott transi- 
tion in a clean system. Moreover, we found that the dis- 
connected term of the density-density correlation func- 
tion, i.e., linig-^o ( n -q)( n q), acts as an easily accessible 
order parameter for the Anderson-Mott transition. We 
found that electron-electron repulsion partially screens 
disorder: for strong interaction electrons feel an effective 
weak disorder potential that should imply an interaction- 
increased localization length. However, once the interac- 
tion exceeds a critical value, a gap opens and the model 
turns into a Mott insulator. From our numerical calcula- 
tions, the ground state is always insulating, yet, upon 
increasing the strength of interaction, the localization 



length may have a non-monotonous behavior when we 
consider the full variational wave function. 

When magnetism is allowed, a compressible and mag- 
netic Anderson insulating phase appears between the 
compressible paramagnetic Anderson insulator and the 
incompressible magnetic Mott insulator. When mag- 
netism is not frustrated, all transitions are likely to be 
continuous. On the contrary, when frustration is included 
by means of next-nearest-neighbor hopping, the param- 
agnetic to magnetic transition turns first order. More- 
over, in the magnetic region, it is also possible to stabi- 
lize many paramagnetic solutions with very low energy, 
suggesting a glassy behavior at finite temperature. In- 
deed, all these paramagnetic states have local moments, 
i.e., magnetic sites that would contribute with a finite 
— ks In 2 term to the entropy at finite temperature. The 
fact that, in this simple two-dimensional model, we find 
local moments in the paramagnetic phase may suggest 
that this is a general feature of disordered systems close 
to a Mott transition. 
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